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The introduction of parallel computers has motivated the search for algorithms which run well on 
them. In this search many new "multigrid- like" algorithms have been proposed. Being inherently 
more parallel than standard multigrid, these algorithms have the potential of being efficient on 
massive parallel machines. This paper considers the parallel mnltigrid algorithm of Frederickson 
and McBryan [13]. This algorithm uses multiple coarse grid problems (instead of one) in the hope 
of accelerating convergence. In this paper, we analyze the convergence properties of this new 
algorithm. This analysis reveals a close relationship with traditional muldgrid methods. 
Specifically, the parallel coarse grid correction operator is identical to a traditional muldgrid 
coarse grid correction operator except that the mixing of high and low frequencies caused by 
aliasing error is removed. We show how to choose appr op ri ate relaxation operators to *airt» 
advantage of this property. Comparisons between the standard muldgrid and the new method are 
made. 
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1. Introduction. The multigrid algorithm is a fast, efficient method for solving a 
wide class of partial differential equations and is now used in many areas of scientific 
computing (such as computational fluid dynamics and structural mechanics) [3], [5], [16], 
[19]. Despite advances in both numerical algorithms and computer hardware, many appli- 
cations require still greater performance than is offered by traditional computers. Given 
the success of the serial algorithm, it is natural to consider parallel multigrid algorithms. 

Many successful parallel multigrid algorithms based on domain decomposition have 
been proposed and/or implemented [1], [2], [4], [7], [6], [9], [10], [14], [18], [20], [21], [23], 
[24], [25], [26], [27], [28], [29], [30], [11]. In these algorithms, parallelism is obtained by 
subdividing the physical region and assigning the subdivisions to different processors. The 
parallel algorithm is then identical to the serial algorithm, with each processor updating 
only grid points in its subdomain. Unfortunately when implemented on a massively par- 
allel computer, this algorithm results in many inactive processors. Consider for example, 
a problem partitioned so that each processor contains only one point on the finest grid. 
If an h to 2h coarsening is used on a two dimensional problem, 3/4 of the processors are 
idle when processing the next coarsest grid. Even more processors are idle on still coarser 
levels. This idle processor problem has spurred research into multigrid-like methods suited 
for massively parallel systems. See [8], [12], [13], [15], [17], [22]. 

In this paper we consider a promising new parallel algorithm which avoids the idle 
processor problem, the Frederickson-McBryan parallel multigrid algorithm [13]. This al- 
gorithm (to be described in section 2) uses multiple coarse grid corrections to improve the 
convergence rate over that of the standard multigrid algorithm. The number of coarse 
grid corrections is matched so that the same number of grid points are processed when 
computing on a fine or coarse grid. For the one dimension Poisson equation it is a direct 
method and is in fact equivalent to odd-even reduction. In this paper we analyze the new 
algorithm for a model anisotropic problem. This analysis reveals a close relationship with 
traditional multigrid methods. Specifically, the parallel coarse grid correction operator is 
identical to a traditional multigrid coarse grid correction operator except that the mixing 
of high and low frequencies caused by aliasing error is removed. To take advantage of 
the aliasing removal, we define a new smoothing number (similar to the Brandt smooth- 
ing number) and a corresponding method for determining relaxation operators for the 
method. Comparisons between the new method and the standard algorithm are made. 
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2. Basic Algorithm* For simplicity, we consider only the situation of one grid point 
per processor. The principle idea in the Frederickson-McBryan algorithm is to use other- 
wise idle processors to perform additional coarse grid corrections hnd hopefully improve 
convergence. For example, on a one-dimensional problem using a two-level method, two 
coarse grid corrections are performed in parallel. Half the processors produce one cor- 
rection by projecting the fine grid equations on the odd points and the other processors 
produce another correction by projection on the even points. The hope is that by com- 
bining these coarse grid corrections the convergence rate can be improved. 

Below we summarize a ‘V’ cycle variant of a 2-dimensional algorithm using one post- 
relaxation sweep to solve the problem 

(1) A\u = 6. 


proc multi( A,,tx,6, level ) 

{ 

if ( level = CoarsestLevel ) then ti = Af b 
else 

Relax(6, ujevel) 

ComputeResidual( Ai, u, b, level, res) 

Pro jResidual(res ,Too , Tc^t go, Teg, level) 
multi( Aooy ^,^ 00 , level+1) 
multi( Ace, *<*» t*oe> level+l) 
multi( Aeo, reo, Ugo, level+l) 
multi( Age , fee , tiee , level+ 1) 

Interpolate(u 00 ,u oe , Ugc,Ugg, level, correction) 
u = u+correction 
endif 

} 

Ai is the current fine grid operator, and Af is the corresponding pseudo-inverse of the 
fine grid operator. 

Age is the coarse grid operator defined on the even points of the current fine grid, 
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Aoo is the operator defined on the odd points of the current fine grid, 

j4oe is the operator defined on the odd points in the x direction and even points in the y 

direction, and 

Aeo is the operator defined on the even points in the x direction and odd points in the 
y direction. It is further assumed that the same projection and interpolation stencils are 
used for each correction (distinguishing this algorithm from the highly parallel multigrid 
algorithms of Hackbusch [17] and Douglas- Miranker [12]). 

The propagation of the error, e*, for a 2 level method can be described by the matrix 


T: 


(2) 

e k = Te fc _x 

where 


( 3 ) 

T = 5(7 - C), 

( 4 ) 

5 = p(G) = l + £>(7, 

i=i 

and 


( 5 ) 

C = PA^RAx. 


S represents the smoothing operator and is a matrix polynomial with coefficients p,. 
For our discussion the iteration operator, G is simply A\. 

C represents the coarse grid correction operator. 

P is the composite prolongation (interpolation) operator (applied at all points). 

R is the composite restriction (projection) operator (applied at all points). 

Ai is the composite coarse grid operator applied at all points (i.e. a combination of 
and Ago ). 

A\ is the fine grid operator. 

The matrices defined above are square ( n 2 x n 2 entries for an n x n grid ) as each 
operator is applied to every grid point. This is in direct contrast with the standard 
algorithm where some operators are applied only to a subset of the points resulting in 
rectangular matrices. 
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The principle difficultly with the above algorithm is to determine smoothing, interpo- 
lation and projection operators to produce convergence rates that are significantly better 
than the standard multigrid approach. The remainder of the paper will address this topic. 

3. Fourier Analysis Framework. The model problem analysis of the parallel al- 
gorithm is somewhat easier than that of standard multigrid. This is primarily because 
the Fourier transform diagonalizes the iteration operator. We illustrate the approach for 
a two dimensional problem defined on the unit square with periodic boundary conditions: 

( 6 ) U xx + flUyy = /. 

The basic analysis for the Poisson equation (/? = 1) appears in [13]. Discretization by 
central differences for an n x n grid yields the linear system: 

(7) Am = h*f 

with h = 1/n. From the previous section, the iteration operator is given by: 

(8) e k = p(Ai)(I - PA%RAi)e k -\. 

A simplification occurs if we assume that the matrices P and R are circulant (i.e. the 
same stencil is applied at each point in the domain). Since Af and A\ are also circulant, 
all the matrices in (8) co mm ute. This implies that the individual choices of P and R are 
unimportant, only the quantity: 

(9) Q = PR. 

With this simplification the error is given by: 

(10) e k = p(A 1 )(I-QAtA 1 )e k - l 


For the remainder of this paper we assume that the restriction operator is simply 
point-wise injection (i.e. R = I), and the interpolation operator (hence Q) is given by a 
local symmetric 9 point stencil. 


( 11 ) 


Q = 


9n 9i 9u > 
9i 90 9i 
\ 9n 9i 9n ' 
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We do not consider stencils of greater width, or whether larger stencils yield significantly 
better convergence rates. 

We transform the operators into Fourier space by considering 

(12) «(*,!/)= u k je 2m{kx+Jv) where P=t;- 

i<IM<p 

In Fourier space all the operators axe diagonal with elements: 



A i 

= 2(1- 

- Xi) + 2f3(l — x 2 ), 

(13) 

A 2 

= (1- 

*i) + 0(1 - xl), 


Q 

= 9o + 

29 i(xi + X2) + 4511X1X2. 

where 




(14) 

Xl = 

cos(27rfc/n) and X2 = cos(2jrj/n ) 


Notice that xi < 0 (xi > 0) corresponds to high (low) frequency in the x direction and 
that X 2 < 0 (®2 > 0) corresponds to high (low) frequency in the y direction. For now we 
take the coefficients p, in the polynomial p(x) and g,- in Q as unknowns. Using these we 
can write down an explicit expression for the error in the parallel algorithm and obtain 
an ‘optimal 1 method by choosing the free parameters to minimize the error over the entire 
range. That is 

(15) 
where 

(16) f(xi,x 2 ) = p(Ai)[I-QAl[Ai]. 

In [13], IVederickson and McBryan used error expressions derived for their parallel multi- 
grid algorithm combined with an optimization routine to deduce optimal parameters. 
With these parameters, they illustrate that extremely fast convergence can be obtained. 
In the next sections the analysis is extended to deduce near optimal parameters analyti- 
cally. Before proceeding with the smoothing, we explain why the parallel method yields 
faster convergence rates than the serial one. 


min 

gj# — 1< X 1 <1 


max _ |T(xi,X2)| 
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4. Cancellation of Aliasing with Multiple Coarse Grids. In this section, we 
show that it is possible to eliminate most of the aliasing error associated with standard 
multigrid methods by using multiple coarse grid corrections. The result is shown for a one 
dimensional problem for simplicity. The same arguments extend to higher dimensional 
problems. 

We first state a le mm a, which is used in the following theorem. 


Lemma 1. 

71 

Ro(avi + 0v h ) = avi - /?tL/ for |/| < 

Rciavt + 0v h ) = 0 ^ + fiv-t for |/| < 

Rl{avi + 0v- /) = a{vi - v - h ) + 0{v-i - v h ) for |/| < 

Re(avi + 0i Li) = a(v { + v_/,) + + v h ) for |/| < 


where 





[vfclj = e 2wik{i/n) 

i f l W 

1*1 = *2 

and j = 1, 

(18) 

= e 4 ” fc(i/n) 

1*1 = 1,..., £ 

and j = 1, 


[{»*], = € 2«(2i-l)(>/r.) 

1*1 = 1— ^ 

and j = 1, 


Re and Rg denote injection onto the even and odd points respectively, and h =% — l when 
l > 0 and h = ^ + / when l < 0. 

Proof 1. For l > 0: 


(19) 


and 


( 20 ) 


Re{av, + 0v h ) = ae 4 '*'W"> + 0e 4irih ^ 

= at* + /?e 4 ” J < n/2_,)/ " 

= otSl + 

= otvi + /3v_i. 

R 0 (av t + 0v h ) = ae 2wilw ~ l),n + 0 e 2wih{2i ~ l),n 
= avi + /3e x ’( n-2, )( 2j-1 )/ n 
= at)/ + /3[c -2 '*^ 2j ~ 1 )^ n c^ 2,_1 ^] 
= a»i — /?v_j 
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We summarize R? and Ro by 


( 21 ) 


( 22 ) 


where 


R e = (v, 


Ro = (Vi 


V-l) 




0 

-I 


0 

I 


0 

I 



(23) 


vi € V t , Vh € K, v-i 6 V_j, 6 V. h , 


(24) vi € Vj, v-i G VL/, t>i G Vj, and i-i € V-/. 

The properties of Rj and Rf follow from (21) and (22). 


The lemma states how the I th Fourier mode on the fine grid (t?j), transforms when it is 
projected onto the coarse grid. Specifically, if the I th mode corresponds to a low frequency 
(|/| < n/4) , then it appears as the I th Fourier mode on the coarse grid. However, if the 
mode ( vh ) corresponds to high frequency (\h\ > n/A ), then it is aliased and appears as 
the —I th mode on the coarse grid. The essential point is that the aliasing on even and odd 
grid points is of opposite sign. 

We illustrate the cancellation of aliasing by defining two coarse grid correction meth- 
ods: single coarse grid correction (SC) tfhich solves on the even points, and multiple 
coarse grid correction (MC) which solves on both the even and odd points as described in 
section 2. 

(25) Csc = P\AfR\ 
and 

(26) C MC = .5 * (PiA+Ri + P 2 A+R 2 ) 
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where Pi , R\ , P2> ^2 correspond to interpolation and restriction on the even and odd points 
respectively. Ae and Ao denote the coarse grid differential operator on the even and odd 
points respectively. We assume that the interpolation and restriction operators employ 
the same stencil for both coarse grids. Notice that the multiple coarse grid correction 

(26) is simply the average of two standard coarse grid corrections on both the even and 
odd points. For the theorem that follows, we assume that the operators Ae and A 0 are 
identical except shifted to operate on different points. Additionally, we assume that the 
operators treat a positive (/ > 0) and negative (/ < 0) mode identically. These assumptions 
corresponds to the following general representations of the pseudo-inverses for arbitrary 
Vi and in 

n n 

(27) A+vi = and A+vj = 

.=1 i=l 

with the same coefficients for both operators. We note that all these assumptions 
are not unre alis tic and hold for periodic constant coefficient operators. 

Without loss of generality, we write the interpolation and restriction operators as: 

(28) Pi = ZRj, P 2 = ZRl, Ri = ReY and R 2 = R*Y 

where Re and Rq are injection operators on the even and odd paints, and Z and Y are 
general n x n circnlant matrices. A property of drcnlant matrices is that the V/’s defined 
in the lemma are eigenvectors. We denote the corresponding eigenvalues as zi and y;. 
That is 

(29) Zvi = zivi and Yvi = yivi. 

We now state the theorem. 

THEOREM 1. Consider a splitting of the Fourier modes (V ) into high and low fre- 
quency: 

(30) V = [V h V h ] 
where 

(31) V t = span 1 <| fc |<„ /4 {vj k ), and V h = spon n/4 <| fc |< n/2 {t; J t}. 
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Using this Fourier mode splitting, we define a portioning of the operator Csc in Fourier 
space: 



where 


(33) C S c = [V l ,V h ]Csc[Vi,V h ] T . 


Then, Cmc is given by: 



where 

(35) C MC = [Vi, V h )C M c[Vi,V h } T . 


Proof 2. 

We look at a matched pair (l, h) from the low and high frequency space: 

(36) vi 6 Vj and v/, G V^. 

Then 

C, sc(<*vi + fivh) = PiA+ Ri(avi + fivh) 

= ZR?AtReY(a Vl + p Vh ) 

(37) = ZR*A+R e (ay t v l + Py h v h ) 

= ZR'A+(ayivt + 0y h v-i) 

n 

= ZRj ( ayivi + fiyhV-i) £ *<1,1. 

*=1 


Csc(avi + fivh) = Z(ayi(vi + v-h) + Pyh(v-i + vh)) Y'. fc.w 

«= 1 

” n 71 

(38) = ayi ^2 k m 2 i v i + PVh Y1 k i\i\ z -i v -i + a Vi 5Z k m z -hv-h+ 

«=1 »=1 1=1 

n 

Pvh Y1 k m z ^- 

*=1 


10 



A similar procedure for the odd points yields: 


(39) 


P 2 A+R 2 (av t + 0v h ) = ayi^kni\zivi - /3y h - ay, ^ ki\,\z_ h v_ h + 


*=1 


1=1 


1=1 


&Vh S *t|/|*A»h- 


i=l 


Averaging ($5) and (SP) yields: 


n n 

(4°) C MC (av { + /3v h ) = ayt £ 

.=1 »=i 

Notice that the av_/, and the /3v~t terms in (88) and (89) cancel when they are av- 
eraged. Additionally, the first term in (88) corresponds to F\, the second term to F3, the 
third term to F4 and the fourth term to F 2 . It is clear that the multiple grid method yields 
the identical expression with the terms from F3 and F4 removed. 


The theorem states that for every multiple grid method of the type described in section 
2, there is a closely related single grid method. In particular, the multiple coarsening 
correction (MC) is identical to the standard one (SC) except that the mixing of low 
and high frequencies (F3 and F4) is removed. These terms F3 and F4 correspond to an 
artificial mixing of high and low frequencies introduced by the method. Typically, these 
terms degrade the performance of the multigrid algorithm as the aliased frequencies on 
the coarse grid bear no physical relationship with the original frequencies on the fine 
grid. While F3 and F4 are error terms, F\ is responsible for the fast convergence of 
the multigrid algorithm. That is, the coarse grid correction accurately reflects the low 
frequency behavior of the fine grid. Since the multiple coarsening method retains the low 
frequency behavior of the single grid method without the aliasing error, it seems logical 
that the multiple coarsening method will converge faster than the standard method. 

5. Interpolation. In this section, we motivate our use of a scaled bilinear interpo- 
lation operator for the new method based on the close connection between the new and 
standard coarse grid corrections. 

The essential idea is to choose an interpolation operator Q that ensures a fast con- 
vergence rate for the multiple coarsening algorithm. Notice that the standard multigrid 
convergence depends on the low to low frequency operator F\ which in turn depends on 
the choice of Q. Since the low-to-low frequency operator in the MC method is the same as 
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that in the standard method, it seems intuitive to choose the same interpolation operator 
for both methods. Therefore, for the MC algorithm, we are led to the choice of bilinear 
interpolation which has been shown to be a good interpolation operator in standard multi- 
grid algorithms for the solution of second order equations. This leads to the following for 

Q- 


(41) 



.25 

.5 

.25 

.5 

1.0 

.5 

.25 

.5 

.25 


Note that the scaling factor 1/4 is needed because we are averaging 4 interpolants from 
multiple coarse grids. This is in fact the same operator that is used in [13] for the 
Poisson equation and is proved to be the ‘best’ using Fourier Analysis arguments. For the 
remainder of this paper we assume that the interpolation operator is given by (41). Of 
course, the possibility exists that it may be better to use interpolation with wider stencils 
in the new parallel method but we shall not explore that here. 

6. Smoothing Criteria/Richardson Relaxation. As we shall see, the typical 
relaxation criteria for standard multigrid methods is not entirely appropriate for the 
Frederickson-McBryan method. In this section we develop a new criteria for the parallel 
method by making use of the close relationship between the standard and new methods. 
We begin by studying Richardson relaxation as a smoothing operator within multigrid 
methods. 

A general Richardson relaxation method can be described by a matrix polynomial. 
For example, the Richardson method for the model problem in section 3 can be written 
as: 

n 

(42) S = P(Ai)= 1 + 

i= 1 

Given such a scheme, the key question is how to determine the p^s so that the multigrid 
method converges rapidly. Probably, the most popular criteria for these parameters is 
that they minimize the Brandt smoothing number [3]. For our model problem, the Brandt 
number is defined as follows: 


(43) 


Hb = max |5(ii, 12)! for *1 < 0 or *2 < 0, 
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where 5 in Fourier space is: 


fl 

(44) 5(ii,i 2 ) = 1 + ^P.[Ai(*i»*2)]' 

«=i 

and xi,x 2 , and A\ are defined in section 3. The Brandt number measures the least 
that any high frequency error is damped after the relaxation is applied. The intuition is 
that the coarse grid correction effectively solves in the low frequency range. Thus, the 
smoothing operator must only damp the error in the high frequency range. Smoothing 
parameters deduced by minimising the Brandt smoothing number are usually easy to 
compute for model problems and work well in multigrid algorithms. In fact, for many 
problems (inclu ding our model problem) the convergence rate for a standard two-level 
multigrid method (for example using point-wise injection and bilinear interpolation) is 
approximately equal to the Brandt smoothing number and is close to optimal over all 
possible p/s. 

Example 1 . Consider the model problem 


(45) 


U xx + /3tiyy — / 


with periodic boundary conditions defined on the unit square. The Brandt number for this 
problem can be written as: 

n 

(46) fi b = max |5(Ai)| = |1 + ^p,Ai| for 2 < A\ < 4 + 40. 

t=i 


where the Brandt number ranges given in ( 4 8) have been transformed into functions of A \ . 
The minimization of (46) over all pi ’s is given by a Tchebyckeff polynomial. For n — 1, 
the optimum is: 

(47) r( i l)=1 -_iL_ 

with 


(48) 


-1 

91 ~ 3 + 2/3 


and fib = 


2/9 + 1 
3 + 2)8' 


It is natural to ask how well this choice of pi performs for a multiple correction 
method. Specifically, consider a two-level single correction method using one iteration of 
this smoothing criteria, point-wise injection, and bilinear interpolation. We can compare 
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•t. <1.1 «. •.» i. 


FlG. 1 . Dark region indicates areas where weight function, N , is equal to 1 in Fourier space. Note: high 
frequency in both directions corresponds to lower left comer of diagram. 


this with a multiple correction method, whose coarse grid correction is simple the average 
of four single corrections. We omit the detailed results and simply state that there is 
almost no improvement using the parallel method over the serial method when the Brandt 
smoothing number criteria is used (i.e. the convergence rates are identical). This is 
because the Brandt smoothing number criteria equi-dis tributes the errors in F 2 ,F3, F4. 
Now that F3 and F4 are zero, we should chose the relaxation operator to minimize the 
effect of 1*2 only. 

To deduce a new smoothing criteria, we examine the motivation behind the Brandt 
number. Specifically, we view the Brandt criteria as a simplification of the general op- 
timization problem. That is, we wish to minimize the spectral radius of the multigrid 
iteration operator. The optimal parameters, p,’s, are given by 


( 49 ) 


minp(T) = minp(p(Ai)(I - C)). 
p|* 


One can view the coarse grid correction, J — C, as a somewhat complex weight function 
applied to p(Ai). We therefore, replace it by a simplier weight function, N. In Fourier 
space this operator is defined as follows: 


( 50 ) 


N(x 1 ,x 2 ) 


1 


0 


if ii < 0 or *2 < 0 
if x\ > 0 and ®2 > 0 
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Fig. 2. Damping of Fourier components with the multiple coarse grid correction. High frequency in both 
directions corresponds to leftmost comer in the diagram where the damping is equal to 1 . 

and is depicted in figure 1. The optimization problem with N replacing C is equivalent 
to minimizing the Brandt smoothing number. The definition of N is motivated by the 
behavior of smoothing operators as well as the behavior of coarse grid corrections. Specif- 
ically, the coarse grid correction provides no improvement over the high frequencies and 
so this error must be entirely damped by the smoothing operator. Second, almost all 
damping of low frequencies is effectively handled by the coarse grid correction. 

Unfortunately for the multiple grid method, the Brandt smoothing number is not 
entirely appropriate. However, we can define a new smoothing number by simply choosing 
a new weight function based on the behavior of the multiple grid correction. Intuitively, 
the removal of the aliasing error should significantly improve the coarse grid correction 
over the high and middle frequencies because they contribute more to aliasing error than 
the low frequencies. Smoothing is still needed to compensate for the fact that the coarse 
grid operators do not accurately reflect the fine grid operator over the high frequencies. 
However, the removal of the aliasing error implies that the smoother should focus on 
damping the highest frequency modes (which are most poorly represented on the coarse 
grid). Specifically for our model Poisson problem, the damping of the Fourier modes with 
the coarse grid correction (depicted in figure 2) is given by: 

(51) A (J - C(l,,*j)) = 1 - i(l + U)(l + 
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X. 


FlO. 3. Dark region indicates areas where weight function, N, is equal to 1 in Fourier space. High frequency 
in both directions corresponds to lower left comer of diagram. 

where A() denotes the damping of a Fourier mode due to the coarse grid correction. 
Further, we observe that 

(52) A (I - C(xi,Z 2 )) = 1 when ii = -1 or x 2 = -1 
and 

(53) X(I - C(x\,X 2 )) < 1 for xi > -1 and x 2 > — 1. 

That is, the coarse grid correction even damps most of the high frequency errors. This 
is a consequence of the absence of aliasing error (which degrades the performance of the 
coarse grid correction). However, since the coarse grid correction does not damp the 
highest frequency modes (xj = —1 or i 2 = — 1), all damping of these components must 
come from the smoothing operator. Based on these observation, we define a new weight 
function corresponding to figure 3 given by: 

r 

1 if xi = — 1 or X 2 = —1 

(54) N(xi,x 2 )=l 

0 if xi ^ 1 and X 2 ^ —1 

k, 

and a corresponding parallel smoothing number: 

p p = maxS(xi,i 2 ) xj = — 1 or x 2 = — 1. 


(55) 
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Assuming that the coarse grid correction damps low frequencies sufficiently, we expect 
that the multigrid convergence rate will be close to the smoothing rate. We emphasize 
that this new smoothing number is a heuristic based on numerical experimentation. 

We demonstrate the use of this new smoothing number with a few examples. 
Example 2. n = 1 (damped Jacobi smoothing). 

We determine the extremal values of A\ over the intervals given by x\ = — 1 and 
®2 = —1. For our model problem, the parallel smoothing number is then given by: 

(56) n P = max |l + piAi| 

4<Ai<4+40 

which is minimized by the Tchebycheff polynomial 

(57) T(A 1 ) = *[A 1 -(4 + 2/3)] 


with 

(58) J >1 = * = and ^ = 2 + r 

Example 3. n = 2 (two step smoothing). 

The parallel smoothing number is given by 

(59) fi p = max |1 + PiAi + P 2 ^?|. 

4<Ai<4+40 

which is minimized by the Tchebycheff polynomial: 

(60) k[A\ - (8 + 4/3)Ai + 2/3 2 + 16/3 + 16). 


Therefore 

(61) 

and 


P2 = k = 


1 

2/3 2 + 16/3 + 16’ 


- 8 + - 4/3 

^ 2/3 2 + 16/3 + 16’ 


P 2 

/3 2 + 8/3 + 8 * 

Tables 1 and 2 illustrates the results for the Frederickson-McBryan method using the 
smoothers corresponding to these examples. Specifically, the tables compare the two- 
level parallel multigrid convergence rate using this smoothing with one obtained using a 
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p pi 

Mp( MC ) 

p(MC) 

optimum p{ MC ) 

1. -.16667 

.3333 

.3333 

.3333 

2. -.125 

.5 

.5 

.498 

3. -.1 

.6 

.6 

.598 

4. -.08333 

.6667 

.667 

.664 


Table 1 


Comparison of convergence rate of parallel method which minimizes parallel smoothing number vs. the 
optimized parallel method (n = 1) on a 64 x 64 grid. 


\L 

?! 

P2 

M P ( MC ) 

P{ MC) 

optimum p( MC ) 

p 

-.352941 

.0294118 

.058824 

.103 

.0938 

B 

-.285714 

.0178571 

.142857 

.143 

.142 

fl 

-.243902 

.0121951 

.219512 

.220 

.217 

B 

-.214286 

.0089286 

.285714 

.286 

.283 


Table 2 


Comparison of convergence rate of parallel method which minimizes parallel smoothing number vs. the 
optimized parallel method (n *= 2) on a 64 x 64 grid. 


numeric optimization routine which chooses the relaxation parameters, to minimize 
the overall convergence rate of the two level multigrid process. 

Not only does the new smoothing number produce near optimal relaxation param- 
eters for the model problem, but it also predicts the convergence rate (similar to the 
Brandt number for the standard multigrid method). With the exception of (3 = 1, the 
convergence rates and the smoothing number are nearly identical. When (3 = 1 however, 
the assumption that the coarse grid correction sufficiently damps the low frequencies is 
not valid. That is, the smoothing damps the high frequencies better than the coarse grid 
correction is removing the low frequency components. Nonetheless, the convergence rates 
obtained are not significantly worse than the optimal. It should be noted that this same 
phenomena occurs with the Brandt smoothing number. That is, the Brandt number does 
not accurately reflect the convergence rate when the smoothing number is very small. 
Typically when this happens, it implies that more smoothing was done than was neces- 
sary. In other words, the relaxation operator is smoothing the high frequencies better 
than the coarse grid correction is removing low frequency errors. 
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n 

= 1 

n 

= 2 

n 

= 3 

p 

p( SC) 

p(MC) 

p( SC) 

>>(MC) 

p{ SC) 

p(MC) 

B 

.6 

.333 

.220 

.103 

.0739 

.074 


.714 

.5 

.342 

.143 

.148 

.068 


.778 

.6 

.434 

.220 

.215 

.0740 

B 

.818 

.667 

.503 

.286 

.275 

.111 

5. 

.846 

.714 

.558 

.342 

.327 

.148 

6. 

.867 

.75 

.601 

.392 

.373 

.182 

7. 

.882 

.778 

.637 

.434 

.413 

.215 

8. 

.895 

.8 

.667 

.471 

.448 

.246 

9. 

.905 

.818 

.693 

.503 

.480 

.275 

10. 

.913 

.833 

.715 

.532 

.508 

.302 


Table 3 


Comparison of convergence rates for serial ( with optimum Brandt smoothing number ) and parallel multigrtd 
(with optimum parallel smoothing number ) for a 64 x 64 grid. 


7. Comparisons with Standard Multigrid. In this section we compare the new 
al g orithm with the standard multigrid algorithm on our 2-dimensional model problem for 
a 64 x 64 grid. The standard multigrid method employs a point-wise injection, bilinear 
interpolation, and a smoothing algorithm which minimizes the Brandt smoothing number. 
The parallel alg orithm uses the same grid transfer operators and minimizes the parallel 
smoothing number. 

Table 3 shows the convergence rates using a 1, 2 and 3 step relaxation algorithm. 
From the table, we notice that the multiple grid method always yields better convergence 
rates than the single grid method. Finally in table 4, we compare the serial and the 
parallel methods using smoothing schemes which result in a convergence rate of one half 
for the serial method. This is accomplished by choosing both the smoothing operator and 
the number of smoothing iterations, and varying (3 to obtain a convergence rate of one 
half. One iteration of the smoothing operator corresponds to the polynomial which mini- 
mizes the Brandt smoothing number for the standard method and minimizing the parallel 
smoothing number for the parallel method. The degree of the minimizing polynomial is 
given by the column ‘degree’ and the number of iterations of the smoothing iteration is 
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degree 

sweeps 

P 

/K sc) 

/KMC) 

1 . 

2. 

1.914121 

.5 

.239161 

1 . 

3. 

3.34732 

.5 

.245292 

1 . 

4. 

4.78521 

.5 

.247372 

2. 

1 . 

3.94949 

.5 

.28261 

2. 

2. 

9.63334 

.5 

.272215 

3. 

1 . 

9.71329 

.5 

.294394 

Table 4 



Comparison of convergence rate of parallel method vs, standard method using different smoothing schemes 
which correspond to a serial convergence rate of one half. 


given by ‘sweeps.’ From the table, we can conclude that when the standard method yields 
a convergence rate of one half, the corresponding parallel method yields a convergence 
rate dose to one quarter for this problem. Obviously, additional tests must be performed 
to fully assess the potential of this parallel method. 

8. Conclusions. We have analyzed the Frederickson and McBryan parallel multigrid 
algorithm and have shown that it can produce convergence rates that are significantly 
better than the standard multigrid method. We have shown that the reason for this success 
is that the mixing of high and low frequencies due to aliasing error is removed. In order 
to take advantage of this, however, the relaxation operator must be appropriately chosen. 
To this end, we defined a new smoothing number which reflects the behavior of the new 
multiple grid correction. In general, the relaxation parameters chosen by this smoothing 
criteria are nearly optimal for our model problem, and result in substantial computational 
savings over the standard method. More tests with variable coefficient problems as well 
as more severely anisotropic problems using more sophisticated relaxation operators are 
planned. 
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